/** \file Random.cc
 * Definitions for the Random number generators
 *
 * (c) Copyright Semih Cemiloglu
 * All rights reserved, see COPYRIGHT file for details.
 *
 * $Id: Random.cc 119 2006-06-21 04:57:25Z semihc $
 *
 *
 */

#include "Random.hh"
#include "Global.hh"

using namespace TC;


TC_IDENT("$Id$");


#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
// RNMX should approximate the largest floating value that is less than 1
#define RNMX (1.0-EPS)


double TC::Random(long& idum)
{
  int j;
  long k;
  static TC_THREAD_LOCAL long iy = 0;
  static TC_THREAD_LOCAL long iv[NTAB];
  double temp;

  // Initialize.
  if (idum <= 0 || !iy) {
    // Be sure to prevent idum==0.
    if (-(idum) < 1) 
      idum = 1; 
    else 
      idum = -(idum);
    // Load the shuffle table (after 8 warm-ups).
    for (j=NTAB+7; j>=0; j--) { 
      k = idum / IQ;
      idum = IA*(idum-k*IQ) - IR*k;
      if (idum < 0) 
	idum += IM;
      if(j < NTAB) 
	iv[j] = idum;
    }
    iy = iv[0];
  }

  // Start here when not initializing.
  k = (idum)/IQ; 
  // Compute idum=(IA*idum) % IM without over-flows by Schrage's method. 
  idum = IA*(idum-k*IQ) - IR*k; 
  if (idum < 0) 
    idum += IM;
  // Will be in the range 0..NTAB-1.
  j = iy / NDIV;
  // Output previously stored value and refill the shuffle table.  
  iy = iv[j]; 
  iv[j] = idum;
  // Because users don't expect endpoint values.
  if ((temp=AM*iy) > RNMX) 
    return RNMX; 
  else 
    return temp;
}


